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We study ground states and elementary excitations of a system of bosonic atoms and diatomic 
Feshbach molecules trapped in a one-dimensional optical lattice using exact diagonalization and 
variational Monte Carlo methods. We primarily study the case of an average filling of one boson 
per site. In agreement with bosonization theory, we show that the ground state of the system in 
the thermodynamic limit corresponds to the Pfaffian-like state when the system is tuned towards 
the superfluid-to-Mott insulator quantum phase transition. Our study clarifies the possibility of 
the creation of exotic Pfaffian-like states in realistic one-dimensional systems. We also present 
preliminary evidence that such states support non-Abelian anyonic excitations that have potential 
application for fault-tolerant topological quantum computation. 

PACS numbers: 03.67.-a, 05.30.Pr, 67.85.-d, 73.43.-f 


I. INTRODUCTION 

The possibility of a fault-tolerant topological quantum 
computation [1H6J based upon topological quasiparticles 
that obey non-Abelian statistics (non-Abelian anyons) 
01 motivated much recent interest in the new systems 
that support such quasiparticles. The idea behind the 
topological quantum computation is that non-Abelian 
anyons could be used to encode and manipulate infor¬ 
mation in a way that is resistant to error. Namely, if 
a quantum system has topological degrees of freedom, 
like non-Abelian anyons, then the information contained 
in those degrees of freedom will be protected against er¬ 
rors caused by local interactions with the environment. 
This provides the possibility of using such systems to per¬ 
form fault-tolerant quantum computation without deco- 
lierence. 

Non-Abelian states of matter also present the funda¬ 
mental intellectual challenge of principle and of experi¬ 
mental realization [Toi - fl2| . The understanding of the ori¬ 
gin and properties of non-Abelian phases is far from com¬ 
plete and is at the frontier of current theoretical research. 
The fundamental objectives are the understanding of the 
interplay between topology and quantum mechanics that 
leads to the formation of non-Abelian phases and the in¬ 
vestigation of new models that have non-Abelian quasi¬ 
particles 0,0. 

Non-Abelian anyons first appeared in the context of 
the fractional quantum Hall (FQH) effect 0|, since FQH 
systems are believed to have a series of exotic non- 
Abelian states. Such states, like the Pfafhan state 
013|, which is the exact ground state of quantum 
Hall Hamiltonians with three-body contact interactions, 
have elementary excitations that are non-Abelian anyons. 
Similar states have also been predicted to occur in cold 
atoms 0|, l2lM26l | . superconductors with p-wave pair¬ 
ing symmetry [J|, hybrid systems of superconductors 


with topological insulators and/or semiconductors Uzh 
Hi], and non-Abelian lattice spin models j32|. 

Although non-Abelian states are associated with two- 
dimensional (2D) systems, analogous states can be found 
in certain one-dimensional (ID) models (3^ - l4l| . Ulti¬ 
mately, such ID non-Abelian states must be braided in 
order to compute. This could be achieved by combining 
ID systems into a 2D network as previously proposed in 
the case of Majorana fermions 0|. Understanding how 
to build states that support non-Abelian defects is an 
important building block towards topological computa¬ 
tion. In this paper we refine a previous proposal for a 
Pfaffian-like state that was proposed as an ansatz for the 
ground state of bosonic atoms subject to three-body infi¬ 
nite repulsive interactions and in a ID optical lattice j37j. 
Although such three-body interactions are rare in nature, 
several experimentally realizable methods have been pro¬ 
posed to realize dominant three-body interactions be¬ 
tween bosonic atoms in optical lattices [371, l43l - |47j . In 
particular, three-body interactions can be efficiently sim¬ 
ulated by mixtures of bosonic atoms and molecules under 
conditions that are achievable with current techn ology in 
systems of atoms and molecules in optical lattices |371 . |43| . 

Therefore, the physical system that we consider is 
a collection of bosonic atoms and diatomic Feshbach 
molecules trapped in a ID optical lattice. Under certain 
experimentally achievable conditions, the system can be 
described by an effective Hamiltonian for bosonic atoms 
with two-body and three-body contact interactions (37ij. 
We study the ground states and elementary excitations 
of the system in the limit of infinite repulsive three-body 
interactions and for a range of values of the two-body 
interaction strength. 

The Pfaffian-like ansatz was originally proposed as an 
ansatz for the ground-state wave function of the system 
in the absence of two-body interactions (37]. However, 
our results show that the Pfaffian-like ansatz wave func- 
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tion most closely corresponds to the exact ground-state 
wave function of the system at some finite value of the 
two-body interaction strength. The results also indicate 
that in the thermodynamic limit this value of the inter¬ 
action strength might be close to the value where the 
system undergoes a quantum phase transition from the 
superfluid state to the Mott insulating state, as previ¬ 
ously found within the bosonization approach |48l |. 

Non-Abelian states of matter order their constituent 
particles following a hidden global pattern that is not 
associated with the breaking of any symmetry [ 2 , |49|. 
This leads to a degeneracy that is not based upon simple 
symmetry considerations and is robust against perturba¬ 
tions and interactions with the environment. Topological 
quasiparticles of such systems exhibit an exotic statistical 
behavior. Namely, the interchange of two identical quasi¬ 
particles takes one ground state into another. If two dif¬ 
ferent exchanges are performed consecutively among the 
quasiparticles, the final state of the system will depend 
upon the order in which these exchanges were carried out. 
This ordering dependence is the reason why such states 
and their quasiparticles are called non-Abelian or non- 
commutative. In addition, the quasiparticles of a non- 
Abelian system are neither fermions nor bosons, which 
motivated the name anyons. 

The Pfafhan-likc states that we consider in this paper 
cannot be characterized by any local order parameter and 
exhibit a global hidden order that is associated with the 
organization of bosons in identical indistinguishable clus¬ 
ters. Indistinguishability between the clusters is achieved 
by symmetrization over the subsets of coordinates of each 
cluster. This symmetrization introduces the possibility 
of topological degeneracy in the space of quasiparticles 
and makes these states potential carriers of non-Abelian 
excitations ITS. [50]. 

Our explicit calculations use variationally optimized 
entangled-plaquette states for systems of up to 60 sites. 
These are benchmarked with exact diagonalization (ED) 
studies of systems of up to 14 sites. 

Using the ED method, we first study ground-state 
properties and elementary excitations of the system for 
small system sizes and with periodic boundary condi¬ 
tions. The Pfafhan-like ansatz was originally proposed 
as an ansatz for the ground-state wave function of the 
system in the absence of two-body_interactions (vanish¬ 
ing two-body interaction strength) {37]. However, our ED 
results clearly demonstrate that the Pfaffian-like ansatz 
better approximates the ground state of the system at 
some finite value of the two-body interaction strength. 
This interaction strength increases with increasing sys¬ 
tem size. Also, the overlap of the exact ground-state 
wave function at such a value of the two-body interac¬ 
tion strength and the Pfafhan-like ansatz wave function, 
decreases more gradually with increasing system size in 
the presence of two-body repulsion than it does in the 
absence of the two-body interactions. 

The ED results thus indicate that in the thermody¬ 
namic limit the Pfafhan-like ansatz wave function most 


closely corresponds to the exact ground-state wave func¬ 
tion of the system at some finite value of the two-body 
interaction strength. This might be close to the value 
where the system undergoes a quantum phase transition 
from the superfluid state to the Mott insulating state, as 
previously found within the bosonization approach [481 ]. 
We also present preliminary evidence that these states 
support non-Abelian excitations required for topological 
quantum computation. 

We further study the ground-state properties of the 
system for larger system sizes. Motivated by the re¬ 
cent success of tensor network methods (5jJ] to numer¬ 
ically simulate a variety of strongly correlated models, 
we use the entangled-plaquette-state (EPS) ansatz, also 
called the correlator-product-state (CPS) ansatz, and the 
variational Monte Carlo (VMC) method [52l46ll j. In the 
EPS approach, the lattice is covered with overlapping 
plaquettes and the ground-state wave function is writ¬ 
ten in terms of the plaquette coefficients. Configura¬ 
tional weights are then optimized using a VMC algo¬ 
rithm. Here, the plaquette coefficients that minimize 
the energy are found using the stochastic minimization 
method [60l - [63j . For small system sizes we find that the 
EPS and VMC calculation gives quite accurate estimates 
of the ground-state energy and the one-body and two- 
body correlation functions. 

To examine the proximity of the ground state to the 
Pfaffian-like ansatz for larger system sizes, we calcu¬ 
late the one-body and two-body correlation functions for 
the exact ground-state and the Pfaffian-like ansatz wave 
functions and compare their asymptotic behavior. Since 
the EPS wave function gives quite accurate estimates of 
the correlations within any plaquette, we estimate the 
asymptotic behavior of the correlation functions from the 
values of the correlation functions for the lattice sites 
within a plaquette. 

We study the ground-state properties of the system 
for the system sizes L = 40 and 60 sites. The results ob¬ 
tained within the EPS and VMC approach are consistent 
with the ED results for smaller system sizes and with the 
results for vanishing two-body interaction strength ob¬ 
tained previously using variational matrix product states 
(MPSs). For the system size L = 60 sites, the maximum 
system size that we have considered, the results indi¬ 
cate that at some finite value of the two-body interac¬ 
tion strength U/t = Uc{L) the exact ground-state wave 
function is still very close to the Pfaffian-like ansatz wave 
function. 

The paper is organized as follows. In Sec. II we intro¬ 
duce the effective three-body interacting atomic Hamil¬ 
tonian for a system of bosonic atoms and diatomic Fes- 
hbach molecules trapped in a ID optical lattice. In Sec. 
Ill we review the theory of the Pfaffian-like states in ID. 
In Sec. IV we present ED results for small system sizes 
and with periodic boundary conditions. The results for 
larger system sizes obtained within the EPS and VMC 
approach are presented in Sec. V. In the final section, 
Sec. VI, we draw our conclusions and discuss possible 
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directions for future research. 


II. EFFECTIVE THREE-BODY INTERACTING 
ATOMIC HAMILTONIAN 

We consider a systems of bosonic atoms and diatomic 
Feshbach molecules trapped in a ID optical lattice. The 
system can be described by the Hamiltonian (37j . Ifkil . [65j 

H = H k + H f + Hi, ( 1 ) 

where 

H k = -t a ^(aja i+ i + h.c.) - t m + h.c.), 


H f = y\5m\mj + 


Uaa f f . 9 / f 

—^-d-CL-CLiCLi + -j={m-CLidi 


h.c.)], 


with t = t a , U 2 = U aa - g 2 / 8, and U 3 = U am 7 2 - The ef¬ 
fective Hamiltonian, is the Hamiltonian for a system 
of bosonic atoms in a ID optical lattice with repulsive 
two- and three-body on-site interactions. 

We further assume that U 3 t a : the limit accessi¬ 
ble in typical setups with 87 Rb atoms [37)]. In the limit 
U 3 — > oo, the Hilbert space is projected onto the sub¬ 
space of states with occupation numbers rn = 0,1,2. 
The bosonic operators subject to this condition, that 
is, the condition (ag i ) 3 = 0 , are referred to as three- 
hard-core bosonic operators and satisfy the commuta¬ 
tion relations [ 03 ^, a^-] = Sij( 1 — §( a 3 ,i) 2 ( a 3 ,i) 2 )- Since 

a 3 i\ n i) = (1 — dni,2)V n i + 1 | n i + 1 }) these operators can 
be represented by 3 x 3 matrices of the form 

/° f 0 \ 

03,i = 0 0 V2 - (4) 

\ 0 0 0 / 


and 


Hi = U a 




U n 




The bosonic operators for atoms and molecules are de¬ 
noted a,; and rrii, respectively. The term Hk describes 
the tunneling processes of atoms and molecules. The 
term Hp is the Feshbach resonance term and the term 
Hi describes the on-site atom-molecule and molecule- 
molecule interactions. Here t a , t m , U aa , U am , and U mm 
are hopping matrix elements for atoms and molecules and 
the on-site atom-atom, atom-molecule, and molecule- 
molecule interaction strengths, respectively. The energy 
offset between open and closed channels in the Fesh¬ 
bach resonance model is denoted <5, and g is the coupling 
strength to the closed channel. 

We further assume U aa , U am , Umm > 0, and 8 > 0. In 
the limit y 2 = g 2 /2S 2 <C 1 , the formation of molecules is 
highly suppressed and the effective Hamiltonian for the 
system can be obtained, to first order in 7 , by projec¬ 
tion of the Hamiltonian, <m>, onto the subspace with no 
molecules. The resulting effective Hamiltonian is |37j 


H e ff = -t o y^(qfoi + h.c.) + Ugml 2 y ^( Q !) 3 ( Q ») 3 


- tml 2 [(«i) 2 ( a i+l) 2 + h ' C - 
i 

+ ( Uaa ~ g 2 /S)Y( a i) 2 ( a i) 2 - 


( 2 ) 


In the limit f m 7 2 t a , valid in typical experiments with 
87 Rb (?i 3 |, the effective Hamiltonian, CD. further reduces 
to 


H eff = ~tY( a i a i +1 +/l.C.) + H 2 ^(4) 2 (aj) 2 
i i 

+ E 4 ^( a I) 3 ( ai ) 3 , 
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In terms of these three-hard-core bosons the projected 
effective Hamiltonian is 

Heff = ~t 'y ^ (03 ^ 03^+1 + h.C.) + — 'y ^ ((^ 3 ^) ( a 3 ,i) ) 
i i 

( 5 ) 

with U = 2U 2 - We also note that within the exper¬ 
imental situation that we study, three-body losses are 
strongly suppressed (the binding energy, and therefore 
the released energy, is not larger than the lattice depth) 
[3^ |. In this paper we study the ground states of the 
Hamiltonian, ([5]), for a range of values of U/t and at 
the fixed average filling factor of one boson per lattice 
site. We compare the ground states to the Pfaffian-like 
ground-state ansatz wave function [37 ]. The properties 
of the Pfaffian-like state are reviewed in the following 
section. 


III. PFAFFIAN-LIKE STATES 

The idea of using symmetrized indistinguishable clus¬ 
ter states as ansatz wave functions for non-Abelian ID 
bosonic liquids was proposed by Paredes, Kielmann, and 
Cirac |37l I38l. [50]. The Pfaffian-like ansatz, inspired by 
the form of the ground state for fractional quantum Hall 
bosons subject to a three-body interaction EM! , was 
originally proposed as an ansatz for the ground-state 
wave function of the Hamiltonian, ©, at U = 0 [37]. 

For bosons in the lowest Landau level subject to 
the three-body interaction potential U 3 ^2 i ^_ : j^ k 8 2 {zi — 
Zj)8 2 (zi — Zk), with Zi = Xi + iyi being the complex co¬ 
ordinate in the 2D plane, the exact ground state in the 
limit U 3 —> 00 is the Pfafhan state E! US, S3, E! 

N /2 N /2 

$ 3 oc 5 t4 {J|(4 - z]) 2 - zj) 2 }. (6) 

i<j i<j 


( 3 ) 





4 






FIG. 1. Schematic of the local projector Pi at a lattice site 
i. The operator Pi projects the two identical local degrees 
of freedom onto a new degree of freedom that is symmetric 
under exchange of the two components. Here the operator 
Pi maps the single-site four-dimensional Hilbert space of two 
species of hard-core bosons, j~ and 4- (red and blue spheres) to 
the single-site three-dimensional Hilbert space of three-hard¬ 
core-bosons (green spheres). 


This state is a symmetrized product of two identical 
Laughlin states [66| . 


JV/2 

*2 (?) 

i<j 

with a =1\4~ Since the Laughlin state of each cluster 
is a zero-energy eigenstate of the two-body interaction 
potential — Zj). three particles can never co¬ 

incide in a state of the form of 0. The operator that 
symmetrizes over the two virtual subsets of coordinates 
{zj} and {4} is denoted 

An ansatz for the ground state of the Hamiltonian, (0 , 
at U = 0 was proposed |37]J in direct analogy with the 
wave function, ©: 

N/2 N/2 

^3 oc |sin(xj - a;J)l II l sin (4 - 4)11- ( 8 ) 

i<j i<j 

This ansatz has the same form as the Pfafhan state, ©, 
with the Laughlin state replaced by a Tonks-Girardeau 
state j67j : 


N/2 

^2 oc n \ sin ( x i - x j)\- ( 9 ) 

i<j 

The Tonks-Girardeau state, ©, is the ground state of 
ID lattice hard-core bosons described by the Hamilto¬ 
nian # 2 , 0 - = — t iiA a 2 ,a,i+i + h.c.) and with periodic 

boundary conditions [68]. The hard-core bosonic opera¬ 
tors a 2 ,o-,i obey (a^cri) 2 = 0, allowing only occupation 
numbers of ti? = 0 or 1 boson per site. Here xf = 2-7T /Li 
and i = 1,..., L, with L being the number of lattice sites. 

To write the ansatz wave function in second quantized 
form, we define a projection operator V such that 

l*3>=P(l*£)®l*t>)- (10) 


The projection operator V is a local operator of the form 

V = Vf L , (11) 

where L is the number of lattice sites and V t is the local 
projector at a lattice site i , 

/I 0 0 0 \ 

Vi= oil 0 . (12) 

\ 0 0 0 V 2 ) 

Pi maps the single site four-dimensional Hilbert space of 
two species of hard-core bosons, f an d to the single¬ 
site three-dimensional Hilbert space of three-hard-core 
bosons as illustrated in Fig. 1. 

It can be further shown that the one- and two-body 
correlation functions for the ansatz wave function, m, 
have the asymptotic behavior [37], 68] 

(' a l+A a i ) A~ 1/4 , (13) 

(4+A a i+A a i a i) 

for large A and for a large system size L. The two- 
body correlation function corresponds to the one-particle 
correlation function for on-site pairs [37j . 

( a l+A a l+A a i a *) OC (14) 

(^2l«2, t ,i+A a 2,t,i|^2)(^2l«2,i,i+A a 24,i|^) 

-»■ A- 1 / 2 A- 1 / 2 , 

where (^2 l a 2 , f 7 ,i+A a 2 ,a, l |d , 2 ) A^ 1 / 2 , for large A, is 

the well known result for a Tonks-Girardeau gas [68] ■ In 
other words, although the system of atoms has some kind 
of coherence, with slowly decaying spatial correlations 
oc A -1 / 4 , the underlying system of on-site pairs is in a 
much more disordered state with a fast decay of spatial 
correlations oc A' 1 . 


IV. EXACT DIAGONALIZATION RESULTS 
FOR SMALL SYSTEM SIZES 

We first examine the ground-state properties and ele¬ 
mentary excitations of the system for small system sizes 
and with periodic boundary conditions using the ED 
method. The ground-state properties are calculated for a 
range of values of the two-body interaction strength U/t 
and at the filling factor of one particle per site. We calcu¬ 
late the overlap of the exact ground-state wave function 
of the Hamiltonian, ©, and the ansatz wave function, 
m, the average occupation of sites with one and two 
particles, and the one- and two-body correlation func¬ 
tions. 

The overlap of the exact ground-state wave function 
and the Pfaffian-like ansatz wave function for the system 
sizes L < 14 and as a function of the two-body interaction 
strength U = U/t is shown in Fig. [2] The results clearly 
demonstrate that the Pfafhan-like ansatz wave function, 






5 


m, is a better ansatz for the exact ground-state wave 
function of the Hamiltonian, ©, at some finite value of 
the two-body interaction strength Uc(L) than it is for the 
exact ground-state wave function at U = 0 as suggested 
previously [37]. 

Also, the overlap decreases more gradually with in¬ 
creasing system size L at Uc(L) than it decreases at 
U = 0 (Fig. [3]). This indicates that in the thermo¬ 
dynamic limit (L —> oo), the Pfaffian-like ansatz wave 
function most closely corresponds to the exact ground- 
state wave function of the Hamiltonian, m, at some finite 
value of the two-body interaction strength Uc(L). 

The ED results also indicate that the value of Uc where 
the overlap is maximal increases with increasing system 
size (Fig. |4|. It can also be shown that the value of U 
where the system undergoes a quantum phase transition 
from the superfluid state to the Mott insulating state, 
Usf-mi , decreases with increasing system size. The ED 
results thus suggest that the value of Uc approaches the 
value of Usf-mi with an increase in the system size L 
and that the Pfaffian-like state might be the state at 
the superfluid-to-Mott insulator boundary as previously 
found within the bosonization approach [48| . 

We have further calculated the average number of sites 
with one particle, m, and with two particles, 712 , at the 
filling factor v = 1 and for a range of values of the two- 
body interaction strength U = U/t. The average number 
of sites with one particle and with two particles is 

1 L 

n\ = (15) 

^ i= 1 
1 L 

n 2 = — ^(aJaJaiOi), 

^ i =1 

where rii = oJa* and ni + 2n 2 = v = 1. The values 
of i and ri 2 for the Pfaffian-like ansatz wave function, 


1.00 


0.90 - 

0.00 0.05 0.10 0.15 0.20 0.25 

1/L 

FIG. 3. Overlap of the exact ground-state wave function of 
the Hamiltonian, ©, and the ansatz wave function, m , at 
the filling factor v = 1 for the values of the two-body inter¬ 
action strength U/t = 0 (red symbols) and U/t = Uc(L), 
where the overlap \(ipExact\tpAnsatz}\ is maximal (blue sym¬ 
bols). Dashed gray lines are included as guides for the eye 
and were obtained by the extrapolation of the ED results for 
L < 14. 

m, and for the exact ground-state wave function of 
the Hamiltonian, ©, at U/t = Uc{L ) (where the over¬ 
lap | (ipExact\ipAnsatz)\ is maximal) are shown in Fig. [5] 
The values of n\ and n -2 for the exact ground-state wave 
function are very close to the values of ?7i and ri 2 for the 
Pfaffian-like ansatz wave function, as can be clearly seen 
in Fig. [5] 

We have also calculated the one-body and two-body 
correlation functions, 

Ci = (<4 +a ( 16 ) 

C 2 = {al +A a] +A aiai), 




FIG. 2. Overlap of the exact ground-state wave function of 
the Hamiltonian, m, and the ansatz wave function, GOD , as a 
function of the two-body interaction strength U/t for system 
sizes of L < 14 sites. Here the filling factor v = N/L = 1, 
with N being the number of particles. 


FIG. 4. The value of the two-body interaction strength U/t 
where the overlap \(ipExact\ipAnsatz)\ is maximal (Uc) for sys¬ 
tem sizes of L < 14 sites and at filling factor v = 1. Dashed 
gray lines are included as guides for the eye and were obtained 
by the extrapolation of the ED results for L < 14. 
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0.00 0.05 0.10 0.15 0.20 0.25 

1/L 

FIG. 5. Average number of sites with one particle rn = 
T Ee=i(ni(2 — m)) (blue and red symbols) and two particles 
ri 2 = 2 ^ (g reen and purple symbols) for the 

exact ground-state wave function of the Hamiltonian, © , at 
U/t = Uc(L) (blue and green symbols) and the ansatz wave 
function, (HOI) , (red and purple symbols). Here the filling 
factor v = 1 and Uc(L) denotes the value of the two-body 
interaction strength where the overlap \(ipExact\tpAnsatz)\ is 
maximal. Dashed gray lines are included as guides for the 
eye and were obtained by the extrapolation of the ED results 
for L < 14. 



Ci(Ansatz) Ci(Exact) C 2 (Ansatz) C 2 (Exact) 



FIG. 6. The one-body (blue and red symbols) and two-body 
(green and purple symbols) correlation functions for the ex¬ 
act ground-state wave function of the Hamiltonian, © , at 
U/t = Uc{L) (blue and green symbols) and the ansatz wave 
function, HOI) (red and purple symbols), for the system size 
L = 14. Here the filling factor v = 1 and Uc(L) denotes the 
value of the two-body interaction strength where the overlap 
l(V> Exact\ipAnsatz}\ is maximal. The long-distance scaling of 
the correlation functions Ci(A) oc A” Q1 and C 2 (A) oc A 
is a\ ~ 0.194 and a 2 ~ 0.776 for the exact ground state and 
aa « 0.212 and a? 2 ~ 0.854 for the ansatz wave function. 


for the exact ground-state wave function at U/t = Uc{L) 
and for the ansatz wave function for a system of L = 14 



FIG. 7. Overlap Oi = 0 2 = O , (11811 . for two degenerate 
first excited states of the Hamiltonian, ©. at the filling fac¬ 
tor v = 1 for values of the two-body interaction strength 
U/t = 0 (red symbols) and U/t = Uc(L), where the over¬ 
lap \(il>Exact\tpAnsatz)\ for the ground state of the Hamilto¬ 
nian, ©, is maximal (blue symbols). Dashed gray lines are 
included as guides for the eye and were obtained by the ex¬ 
trapolation of the ED results for L < 14. 


sites. The results in Fig. [G] show that the correlation 
functions for the exact and ansatz wave functions show 
very similar asymptotic behavior. 

Elementary excitations further reveal the topological 
nature of the Pfafhan-like state. By construction, the 
Pfaffian-like ansatz, m, has a hidden global order as¬ 
sociated with the organization of particles in two iden¬ 
tical indistinguishable copies of the same state. Conse¬ 
quently, the elementary excitations above the Pfaffian- 
like ground-state ansatz exhibit non-Abelian statistics 
[50| . The argument for this proceeds as follows. The 
elementary excitations can be constructed by creating a 
quasihole in each of the copies and symmetrizing [5CJ. 
Symmetrization leads to a topological degeneracy in the 
subspace of elementary excitations and non-Abelian al¬ 
gebra of exchanges of elementary excitations (quasiholes) 

If the ground state of the system is close to the Pfaffian- 
like ansatz, elementary excitations above it also exhibit 
non-Abelian statistics. To confirm this statement we 
compare excited states of the Hamiltonian, (©), with the 
ansatz wave functions for the excited states [50J 

\4 m ’ n) ) = n\*l {m) )®\* l 2 {n) }), ( it ) 

where are eigenstates for each of two copies. In 

particular, we calculate the overlap of the first excited 
state of the Hamiltonian, ©, and corresponding ansatz 
wave function for the first excited state, 113. 

For a fixed number of particles N (with N either even 
or odd) an excited state with 2 n elementary excitations 
(quasiholes) that are SU(2) 2 (Ising) anyons (similar to 
Majorana fermions) is expected to have topological de¬ 
generacy 2" -1 [40| |. We find that the first excited state 
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of the Hamiltonian, ©> is twofold degenerate and corre¬ 
sponds to a state with four quasiholes. 

The overlap for each of the two degenerate first excited 
states is calculated by considering the total overlap with 
the manifold of degenerate ansatz states, (ED, that cor¬ 
respond to the first excited states. We find four degener¬ 
ate, linearly independent, ansatz states that correspond 
to two degenerate first excited states of the Hamiltonian, 
©• The degeneracy of the exact first excited state is 
two, and not four, since the Hamiltonian, ©, does not 
have particle-hole symmetry. The overlap is then given 
by 


Oi = ^E i('4 1) i4 L) >i 2 = (is) 

where |V’j 1 ^}> with i = 1,2, are two degenerate first ex¬ 
cited states of the Hamiltonian, ©. and with 

k = 1,..., 4, are corresponding degenerate ansatz states. 
We note that the states \<f>^) form an orthonormal basis 
within the degenerate manifold, which leads to expres¬ 
sion (fT51) for the total overlap. We find that 0\ = O 2 for 
all values of the two-body interaction strength U/t that 
we have considered. 

The overlap of the first excited states and correspond¬ 
ing ansatz wave functions is shown in Fig. [7] In agree¬ 
ment with the results for the ground-state wave function, 
the overlap decreases more gradually with increasing sys¬ 
tem size L at Uc(L) than it decreases at U = 0 (Fig. [3D. 
The results also demonstrate that the overlap with the 
ansatz wave function for four quasiholes is very close to 1, 
indicating that the elementary excitations of the system 
are non-Abelian. 


V. VARIATIONAL MONTE CARLO 
CALCULATION 

Studying properties of the system for larger system 
sizes with the ED method is not possible due to the rapid 
increase in the Hilbert-space size with increasing system 
size. Motivated by successes of the tensor network meth¬ 
ods 0 to numerically simulate a variety of strongly cor¬ 
related models, we further study the properties of the sys¬ 
tem for larger system sizes using the entangled-plaquette- 



FIG. 8. Illustration of the six-site ID plaquettes with a five- 
site overlap between the plaquettes used in the calculation of 
the ground-state properties of the system for larger system 
sizes. 


u/t 

Eo,eps Eo,ed 

R(xl0 2 ) 

0.0 

-1.61170 -1.62516 

0.828 

0.5 

-1.49408 -1.50630 

0.811 

1.0 

-1.38258 -1.39293 

0.743 

« 1.18 « 

U c -1.33919 -1.35329 

1.042 

1.5 

-1.28014 -1.28519 

0.393 

2.0 

-1.17054 -1.18326 

1.075 

2.5 

-1.07458 -1.08740 

1.179 

3.0 

-0.98504 -0.997966 

1.295 


TABLE I. Ground-state energy per site (in units of t) for 
the system size L = 14 and with periodic boundary con¬ 
ditions. Here R stands for the relative error with respect 
to the exact ground-state energy and is defined as R = 
( Eo : ed — Eq^eps) / Eq,ed- 

state ansatz optimized using the variational Monte Carlo 
method [52H6l| . Within the EPS approach, also called 
the correlator-product-state approach, the lattice is cov¬ 
ered with overlapping plaquettes and the ground-state 
wave function is written in terms of the plaquette coef¬ 
ficients. Configurational weights can then be optimized 
using a VMC algorithm. 

For a lattice with L sites, an arbitrary quantum many- 
body wave function can be written as 

i^> = E w ni \n\ ,..., tie) — 'y ' H n |n), (19) 

where n = {m,..., til} denotes the vector of occupancies 
and W n is the amplitude or weight of a given configura¬ 
tion n. For the system described by the Hamiltonian © 
rii £ {0,1, 2} for the lattice sites i = 1,..., L. 

In the EPS (CPS) description of a bosonic system, the 
weight W n is expressed as a product of the plaquette 
coefficients over the lattice 

W ni ,..., nL =H C”», (20) 

V 

where n p = {n p 1 ,..., n p {\ is the occupancy vector of the l- 
site plaquette p. In many cases, the qualitative behavior 
of large systems can be described even by plaquettes with 
a small number of sites. The EPS wave function corre¬ 
sponding to the ground-state wave function of the system 
gives reasonable estimates of the ground-state energy and 
short-range correlations. The estimates improve with an 
increase in plaquette size and greater overlap between the 
plaquettes. 

Here we choose six-site ID plaquettes with a five-site 
overlap between the plaquettes as illustrated in Fig. [5] 
The weights W n for the ground-state wave function of 
Hamiltonian © with periodic boundary conditions can 
then be written as 

yyr (J n \ ,ri2 ,...,riQ * ^n2,ri3,...,n7 . . Qn l . .n*, ( 21 ) 

Within the variational MPS approachf37f this plaquette 
choice would correspond to matrices of dimension \ = 
, with L° = 5 being the number of overlapping lattice 
sites. 







Cl (ED) 
Ci (EPS) 
C2(ED) 
C 2 (EPS) 

0.0 0.5 1.0 1.5 

ln(A) 



FIG. 9. The one-body (red and blue symbols) and two- 
body (green and purple symbols) correlation functions ob¬ 
tained using the ED (red and green symbols) and the EPS 
and VMC (blue and purple symbols) methods for the sys¬ 
tem size L = 14 and at the value of the two-body interaction 
strength U/t ~ Uc(L) « 1.18. 


In a VMC algorithm the energy E is written as 

Z a . a 'W:,(n'\H\n)W a 
W#> En |Wn | 2 

= J2 P " E ^ ( 22 ) 


where it is assumed that the wave function | i/j) = 
En^nln) is not normalized, and the local energy E n 
and the probability P n are given by 

^ = ( 23 ) 

n' n 

p l^nl 2 
11 EnW 


The expectation value of any operator O can be expressed 
in the same form by replacing the Hamiltonian H with 
the operator O. The probability P n is never explicitly 
calculated from Eq. (E51) . Instead, for a given set of pla- 
quette coefficients, the energy can be efficiently computed 
using the Metropolis algorithm jH|. 

Within the Metropolis algorithm, used to sample the 
probability distribution, the overall energy can be effi¬ 
ciently computed as an average of the sampled local en¬ 
ergies. In our calculation the total number of atoms N 
is fixed. We start from a randomly chosen initial con¬ 
figuration |n) = \ni,ri 2 , with EiE n i = N and 

rii £ {0,1,2} for i = 1 and then generate via the 

Metropolis algorithm a large set of new configurations by 
replacing n, with m — 1 (if rq > 0) and rij with rij + 1 
(if rij < 2) at two neighboring sites i and j. Starting 
from configuration |n), the acceptance probability of a 
new configuration |n') is given by 


Pa = min 


' IwW 

_|W(n)|2 


1 


(24) 


According to the variational principle, minimization 
of expression (1^1) with respect to the weights gives an 


upper bound of the ground-state energy. The plaquette 
coefficients that minimize the energy can be found by 
using the stochastic minimization method [60H63i . I69I . I70| , 
which requires only the first derivative of the energy with 
respect to the plaquette coefficients, which is given by 


dE 

dCp p 


2 E 



E, 




EL 


(25) 


where the wave function | ip) in Eq. (1221) is approximated 
by the EPS (CPS) wave function and 


A np = 

p 


1 dW n 

w^ac^ 


bp 

> 

Op 


(26) 


with W n given by Eq. m- Here b p denotes the num¬ 
ber of times the plaquette coefficient C p p appears in the 
product, m, for the amplitude W n for configuration |n). 
If the same plaquette coefficient is used for multiple sites 
(e.g., for a translationally invariant choice of plaquettes), 
b p > 1. If each plaquette coefficient is used once, b p = 1. 
Equivalently to the overall energy the first derivative can 
be efficiently calculated using the Metropolis algorithm 
from the same sample used to compute the overall energy. 


U=0 L=40, 60 ffi5»0.232, 0.205 



U=0 L=40, 60 ff 2 «0.797, 0.724 



FIG. 10. The one-body (Ci) and two-body (C 2 ) correlation 
functions for the system sizes L = 40 (red symbols), and 
L = 60 (blue symbols) sites and at the value of the two-body 
interaction strength U/t = 0. Here the asymptotic behavior 
of the correlation functions is Ci = (a\ +A ai) —> A -ai and 
C2 = (ol + A a !+A a * a *) “A “ 2 . 
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L=60 U/t=0, 1.2, 1.3 a-i-0.205, 0.213, 0.218 



ln(A) 

L=60 U/t=0, 1.2, 1.3 ff 2 «0.724, 0.756, 0.764 



FIG. 11. The one-body (Ci) and two-body (C 2 ) correlation 
functions for the system size L = 60 sites and at the values of 
the two-body interaction strength U/t = 0 (red symbols), 1.2 
(purple symbols), and 1.3 (blue symbols). Here the asymp¬ 
totic behavior of the correlation functions is Ci = (aJ +A ai) — ► 
A" ai and C 2 = {a\ +A a\ +A a iai ) A"“ 2 . 

The steps of the VMC algorithm used to calculate the 
ground-state properties of the system are as follows: (i) 
start from the randomly chosen complex values for the 
plaquette coefficients, (ii) evaluate the energy and its gra¬ 
dient vector, (iii) update all plaquette coefficients Cp p 
according to 

9 C”’ - rS(k) • sign , (27) 

and (iv) iterate from (ii) until convergence of the energy 
is reached. Here r is a random number between 0 and 1, 
and S(k) is the step size for a given iteration k. 

In each iteration k, the energy and its derivative are 
estimated from F(k) x L values, where L is the number of 
lattice sites and F(k) is called the number of sweeps per 
sample. In a given sweep each lattice site is visited se¬ 
quentially and a move, n n', to a new configuration is 
proposed by changing the occupancy numbers m —> m — 1 
(for rii > 0) and rij —> rij + 1 (for rij < 2) at two neigh¬ 
boring lattice sites, i and j. Also, to achieve convergence 
and reach the optimal energy value it is important to 
carefully tune the gradient step 5(k). For each iteration 
k, the number of sweeps F is increased linearly, F = F^k, 
and the procedure of evaluating the energy and updat¬ 


ing the coefficients is repeated G = G^k times. The step 
size is gradually reduced per iteration. Here we use a 
geometric form, 5 = SoQ k , with Q = 0.9. 

The number of sweeps per iteration is increased be¬ 
cause the derivatives become smaller as the energy min¬ 
imum is approached and require more sampling in order 
not to be dominated by noise. An increasing G effec¬ 
tively corresponds to a slower cooling rate. Here we take 
Fq = 100, Go = 10, and Q = 0.9. The initial mini¬ 
mization routine is performed with So = 0.5 for 50 itera¬ 
tions. The resulting plaquette coefficients are then used 
as a starting point for a new run of 50 iterations with 
Sq = 0.05. After the minimization is complete the expec¬ 
tation values are calculated by repeating the procedure 
for a single iteration with 0 step size and large F and 
G to obtain more accurate estimates of the expectation 
values. 

It is also important to note that it is more difficult to 
obtain good estimates of the ground-state energies and 
the correlation functions for small system sizes due to 
the presence of the statistical error in the stochastic al¬ 
gorithm. Having a larger number of parameters allows 
the optimization method more freedom in finding the 
minimum energy state and the statistical error can be 
controlled by increasing the system size. 

The results for the ground-state energy for the system 
size L = 14 and for several values of the two-body inter¬ 
action strength U/t are listed in Tab. U and compared 
to the ED results (exact ground-state energy values). As 
can be seen in Table [I] the relative error, defined as 

r= Eq, ed -Eq , eps ) (2g) 

Eo,ed 

does not exceed 1.3% for any value of the two-body inter¬ 
action strength U/t. The EPS and VMC calculation also 
gives quite accurate estimates of the correlation functions 
as demonstrated in Fig. [21 

To examine the proximity of the ground-state wave 
function for larger system sizes to the Pfafhan-like ansatz 
wave function, we further calculate correlation functions 
for system sizes L = 40 and 60 sites. Since the EPS wave 
function gives quite accurate estimates of the correlations 
within any plaquette p, to estimate the asymptotic be¬ 
havior of the correlation functions we calculate the one- 
body and two-body correlation functions, m, for the 
lattice sites within a plaquette p (A = 1, l p — 1). 

The results for the one-body and two-body correlation 
functions at U = 0 and for the system sizes L = 40 and 
60 sites are shown in Fig. [TUI The values of oq and a 2 
that describe the asymptotic behavior of the correlation 
functions, 

Gi = (al +A a i )->A-“ 1 , (29) 

C 2 = ( a ]+ A a l + A aiai ) —> A 2 , 

move away from the values expected for the Pfaffian-like 
state ansatz (oq = 0.25 and a 2 = 1) with increasing 
system size. The results thus suggest that the overlap 
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between the exact ground-state wave function at U = 
0 and the Pfaffian-like ansatz decreases with increasing 
system size. This is consistent with the ED results and 
with the previously obtained results for the system sizes 
L < 40 obtained using the variational MPS |37| . 

The ED results presented in the previous section also 
suggest that the Pfaffian-like ansatz wave function, m, 
better approximates the exact ground state of the Hamil¬ 
tonian, ©, at some finite value of the two-body inter¬ 
action strength U/t = Uc{L) than it approximates the 
exact ground-state wave function at U = 0. The results 
for the one-body and two-body correlation functions for 
the system size L = 60 sites (Fig. [Till show that the val¬ 
ues of a i and a-i increase with an increase in the value 
of the two-body interaction strength U/t = U. In other 
words, at some value Uc{L) the values of a\ and ai will 
be the closest to the values expected for the Pfaffian-like 
ansatz wave function (aq = 0.25 and a 2 = 1). This in¬ 
dicates that the overlap between the exact ground-state 
wave function and the Pfaffian-like ansatz wave function 
is maximal at Uc{L ), in agreement with the ED results 
for smaller system sizes. 

Previous calculations with variational MPS found the 
values of oq and a .2 for the system size L = 20 sites to 
be ai = 0.22 and a 2 = 0.83 for the exact ground-state 
wave function at 17 = 0 and oq = 0.24 and a 2 = 0.99 
for the Pfaffian-like ansatz wave function (37]. The cor¬ 
responding overlap between the exact ground-state and 
the Pfaffian-like state wave functions was found to be 
ss 0.955 [37j. Also, for the system size L = 40 sites at 
U = 0 the overlap is ~ 0.90 [37], which corresponds to 
oq ps 0.232 and «2 ~ 0.797 obtained within our EPS 
and VMC calculation for the exact ground-state wave 
function (Fig. fTOl) and oq s=s 0.25 and ps 1 for the 
Pfaffian-like ansatz wave function. 

It is difficult to determine the exact values of Uc{L) 
from our EPS and VMC calculation. However, for the 
system size L = 60 sites, the maximum system size 
that we have considered, we find that op > 0.218 and 
a 2 ^ 0.764 at U/t = Uc{L). Therefore, based on the re¬ 
sults mentioned in the previous paragraph, we estimate 
that the overlap between the exact ground-state and the 
Pfaffian-like-state wave functions is still good for the sys¬ 
tem size L = 60 and at U/t = Uc{L). 


VI. CONCLUSIONS 

We have studied ground states and elementary excita¬ 
tions of a system of bosonic atoms and diatomic Feshbach 
molecules trapped in a ID optical lattice. Under certain 
conditions, which are experimentally achievable with cur¬ 
rent technology in systems of cold atoms and molecules 
in optical lattices, the system can be described by an 
effective Hamiltonian for bosonic atoms with two- and 
three-body interactions. We have considered the limit 
of infinitely strong three-body interactions for a range of 
values of the two-body interaction strength. The ground- 


state properties of the system were calculated using the 
ED method for small system sizes, and the EPS and 
VMC method for larger system sizes. 

The Pfaffian-like ansatz was originally proposed as an 
ansatz for the ground-state wave function of the effec¬ 
tive Hamiltonian in the absence of two-body interac¬ 
tions. However, our results clearly demonstrated that 
the Pfaffian-like ansatz wave function is a better ansatz 
for the ground-state wave function of the effective Hamil¬ 
tonian at some finite value of the two-body interaction 
strength. This value of the two-body interaction strength 
might be close to the value where the system undergoes 
a quantum phase transition from the superfluid state to 
the Mott insulating state, as previously found within the 
bosonization approach. We also demonstrate that these 
states support non-Abelian excitations required for quan¬ 
tum computation. 

Further work is necessary to find an experimentally re¬ 
alizable model with a ground-state wave function that 
can be even better approximated by the Pfaffian-like 
ansatz wave function. This can possibly be achieved in 
a system with long-range interactions. An additional di¬ 
rection for future research is to consider similar 2D non- 
Abelian models, for example, anisotropic systems consist¬ 
ing of coupled interacting ID wires. Such anisotropic 2D 
lattice models can have interesting non-Abelian Chern 
insulating phases and fractional topological insulating 
phases. 

In order to use non-Abelian states for quantum com¬ 
putation, one must be able to braid them. Although this 
is not possible for a strictly ID system, creating a net¬ 
work of such ID systems connected by T-junctions, as 
suggested previously in the context of Majorana quan¬ 
tum wires [42;], potentially allows this. In the case of 
Majorana quantum wires a T-junction [42| allows for 
adiabatic exchange of two Majorana fermions. Such a 
T-junction has topological and non-topological regions 
that can be controlled by individually tunable gates. In 
principle, similar T-junction networks can be created for 
the bosonic system that we have considered, where topo¬ 
logical and nontopological regions of the network can be 
controlled by tuning the two-body interaction in different 
regions of the T-junction network (for example, by chang¬ 
ing the depth of the optical lattice in certain regions of 
the T-junction). 

In a similar manner to fractional quantum Hall states, 
non-Abelian anyons can be created by creating pairs of 
quasiholes, with one quasihole in each cluster [33, [5C|. 
These non-Abelian anyons are Ising anyons [SU(2)2 
anyons], similar to Majorana fermions, and excitations 
of the v = 5/2 fractional quantum Hall state (Pfaffian 
state). We also note that braiding of SU(2)2 anyons 
alone does not permit universal quantum computation 
Q. However, to obtain a universal set of gates, braid¬ 
ing of SU( 2)2 anyons needs to be strengthened only by a 
single-qubit 7t/8 phase gate and a two-qubit measurement 
p]. Also, Pfaffian-like states obtained by symmetrization 
of two identical copies can be generalized to states ob- 
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tained by symmetrization of k identical copies that sup¬ 
port SU(2)fc anyons (38t [50l] and can be used for univer¬ 
sal quantum computation. For example SU(2)3 anyons 
(like Fibonacci anyons) can be used for universal quan¬ 
tum computation Q. The results presented here thus 
constitute an important step towards understanding gen¬ 
eralized states that support SU(2)fc anyons. 
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